% function estima_out_analysis 

%% 10 In detail Analysis of MODE 
% 10.1 Preliminaries 
% ANA_IN is structure that handles all inputs 
if ~exist('ana_in', 'var'); % ANA_IN does not exist as a variable
    ana_in = []; 
end
[e,ana_in]=ch_field(ana_in,'flag_addz',zeros(1,size(Y, 2)));
ana_in.flag_mixfreq = flag_mixfreq; 
ana_in.flag_gtv     = flag_gtv; 
ana_in.stanames = stanames;
ana_in.shonames = shonames;
ana_in.parnames = parnames;
ana_in.state_sel = state_sel; 
ana_in.modelnames = {'top mode'};
ana_in.ssnames = ssnames;
ana_in.savepath = save_path;
ana_in.flag_acov = 0; 
ana_in.flag_acdecomp = 0; 
ana_in.nper = 3;  
ana_in.vecsim = [100 size(Y,1)]; 
ana_in.flag_medbonly = 0; 
ana_in.nrep          = 4000; 
if flag_mixfreq == 1
    [out_tag,cellrep]=tag_admissable(out.case_tag);
    ana_in.case_tag = out.case_tag; 
else 
    ana_in.case_tag=[]; 
end 
% 10.1 Code for the analysis 
ana_out = lmj_modeanalysis_compute(funcmod, parmode_top, Y, ana_in, solveopt, addsol);
pertmat = ana_out.pertmat;


%%
% =========================================================================
% 11. Volatilties: Simulated Moments and Plot of Histogram 
if flag_mixfreq == 0; 
    stdsimz_mat = ana_out.stdsimz_mat; 
elseif flag_mixfreq == 1; 
    stdsimz_tag_mat = ana_out.stdsimz_tag_mat;
    stdsimz_mat     = ana_out.stdsimz_mat;
end 

%% Create table with the volatilities 
if  ~isempty(stdsimz_mat); % Base Frequency
    for ii = 1:ana_out.num_valid;
        
        %% Cell and Latex Table with output from High Frequency 
        [std_outcell,std_tab]= playmod_outdraws( (squeeze(stdsimz_mat(:, :, ii)))', znames);
        std_outcell = std_outcell(3:end, :); 
        MakeLaTeX_table([std(Y)' std_tab.med std_tab.mean std_tab.p5 std_tab.p95],znames,{'Data','Median','Mean','5th','95th'},'Series,Simulated Moments',['Std, High Frequency T=',num2str(ana_in.vecsim(2)),',',subfname(2:end)],['Tab Std High Frequency ',subfname(2:end)],0,outpath);
        std_data    = ['Data'; num2cprec(std(Y))'];
        std_outcell = [std_outcell(:, 1), std_data, std_outcell(:, 2:end)]; 
        std_outcell = merge_cells({'High Frequency'}, std_outcell); 
        clear std_tab std_data; 
        
        %% Cell and Latex Table with output from Low Frequency
        %  High frequency series will be averaged
        if flag_mixfreq == 1 && ~isempty(stdsimz_tag_mat);
            [std_outcell_tag,std_tab_tag]= playmod_outdraws((squeeze( stdsimz_tag_mat(:, :, ii) ) )', znames);
            std_outcell_tag              = std_outcell_tag(3:end, :);            
            % For the actual data, Case TAG for LF must be 1, since only
            % the end-point is observable, Case Tag 1=end of period, Case tag 2=taking averages                                                                                               
            out_tag_std = tag_admissable([2*ones(out.nzHF, 1); ones(out.nzLF,1)]); 
            Y_LF        = tag_datasim_general(Y, out_tag_std, out.NUMIT); 
            std_LF      = std(Y_LF); 
            std_data = ['Data'; num2cprec(std_LF(:))];            
            MakeLaTeX_table([std_LF(:) std_tab_tag.med std_tab_tag.mean std_tab_tag.p5 std_tab_tag.p95],znames,{'Data','Median','Mean','5th','95th'},'Series,Simulated Moments',['Std, Low Frequency, T=',num2str(size(Y_LF,1)),', ',subfname(2:end)],['Tab Std Low Frequency ',subfname(2:end)],0,outpath);            
            std_outcell_tag =[std_outcell_tag(:, 1), std_data, std_outcell_tag(:, 2:end)];
            std_outcell_tag = merge_cells({'Aggregated Frequency (Quarterly)'}, std_outcell_tag);
            std_outcell     = merge_cells(std_outcell, std_outcell_tag, 1);
        end
        if ii == 1;
            std_outcell_full = std_outcell;
        else
            std_outcell_full = merge_cells(std_outcell_full, std_outcell, 2 );
        end
    end
    std_outcell_full = merge_cells(topcell, std_outcell_full, 1);
    
    cd(outpath)
    if ~isunix
        xlswrite(mode_filename,std_outcell_full,'std_sim');
    else
        savecell(std_outcell_full,[],outpath,['Est Result', '_', mode_filename]);
    end
    cd(cucd)
    
    % 1) plot volatilities
    in_stdhist.modenames = ana_out.modelnames;
    in_stdhist.varnames  = znames;
    in_stdhist.savepath  = ana_in.savepath;
    in_stdhist.savename  = 'Z STD';
    plot_stdhist(stdsimz_mat, Y, in_stdhist);
    
end

%%
% =========================================================================
% 12 Impulse Responses 
% 12.1 For ALL observables 
if ~isempty(ana_out.irfz_mat);
    in_irf.modenames = ana_out.modelnames;
    in_irf.varnames = znames;
    in_irf.shonames = shonames;
    in_irf.savepath = save_path;
    in_irf.savename = 'Z IRF';
    in_irf.addlev = ana_in.flag_addz;
    [e, in_irf] = ch_field(in_irf, 'flag_irfbysho', 1);
    plotmultiple_irfs(ana_out.irfz_mat, in_irf);
end
% 12.2 For selected states 
if ~isempty(ana_out.irfs_mat)
    in_irf.modenames = ana_out.modelnames;
    in_irf.varnames =  ana_in.state_sel;
    in_irf.shonames = shonames;
    in_irf.savepath = save_path;
    in_irf.savename = 'Sel States IRF';
    [e, ana_in] = ch_field(ana_in, 'flag_adds', zeros(size(ana_in.state_sel)));
    in_irf.addlev = ana_in.flag_adds;
    [e, in_irf] = ch_field(in_irf, 'flag_irfbysho', 1);
    plotmultiple_irfs_split(ana_out.irfs_mat, in_irf);
end
% 13) Variance Decompositions.
% 3.1) of the observables
if ~isempty(ana_out.vdechz_mat)
    in_vdec.modenames = ana_out.modelnames;
    in_vdec.varnames =  znames;
    in_vdec.shonames = shonames;
    in_vdec.savepath = save_path;
    in_vdec.savename = 'Z VDEC';
    plot_vdecs(ana_out.vdechz_mat, in_vdec);
end
% 3.2) of the variables
if ~isempty(ana_out.vdechs_mat)
    in_vdec.modenames = ana_out.modelnames;
    in_vdec.varnames =  ana_in.state_sel;
    in_vdec.shonames = shonames;
    in_vdec.savepath = save_path;
    in_vdec.savename = 'Sel State VDEC';
    plot_vdecs(ana_out.vdechs_mat, in_vdec);
end
%
%%


% 4) AUTOCORRELATION
% 4.1) of the obervables;
in_ac.modenames = ['data';  ana_out.modelnames(:)];
in_ac.varnames = znames;
%in_ac.varnames = {'u'; 'Y'; 'C'; 'I';  'v'; 'f'; 's'}; 

in_ac.savepath = save_path;
in_ac.savename = 'Obs AC Median';
if flag_mixfreq == 1
    %=========================================================
    % Now report aggregations for even high frequency data
    % by Xi Luo 04/07/2011
    in_ac.flag_block = 0; 
    
    isHF = out.case_tag == -1; 
    isLF = out.case_tag ~= -1; 
    
    % extracting HF data
    Y_HF = Y(:, isHF); % Original High frequency data
    out_tag_HF = tag_admissable(2*ones(out.nzHF, 1)); % 2 means taking averages
    Y_HF_tag = tag_datasim_general(Y_HF, out_tag_HF, out.NUMIT); % transformed into low frequency
    
    % extracting LF data
    %out_tag_LF = tag_admissable(1*ones(out.nzLF,1)); 
    %Y_LF = tag_datasim_general(Y(:, isLF), out_tag_LF, out.NUMIT); 
   
    
    % constructing data autocovariance
    AC_data_HFtag_LF = nanaccov([Y_HF_tag, Y_LF(:,isLF)], (0:ana_in.nper)); % dealing with NAN
    AC_data_HF = accov(Y_HF , (0:ana_in.nper)); 
    ac_data = AC_data_HFtag_LF;
    ac_data(1:out.nzHF, 1:out.nzHF, :) = AC_data_HF; 
    %=========================================================
    
   
elseif flag_mixfreq == 0
    ac_data=accov(Y,(0:ana_in.nper));
end

if flag_mixfreq == 0; 
    acsimz_mat = ana_out.acsimz_mat;
elseif flag_mixfreq == 1; 
    %acsimz_mat = ana_out.acsimz_tag_mat;
    acsimz_mat = ana_out.acsimz_mix_mat;
end

%% 
ac0data_cell =[['AC(0)', znames(:)'];  [znames(:), num2cprec(ac_data(:,:,1))] ];
%

aczero_cell=out_vdec(squeeze(acsimz_mat(:,:,1,:)),znames,znames,'Model AC(0) Sim'); 
%

ac1data_cell =[['AC(1)', znames(:)'];  [znames(:), num2cprec(ac_data(:,:,2))] ];
%
acone_cell =out_vdec(squeeze(acsimz_mat(:,:,2,:)),znames,znames,'Model AC(1) Sim'); 
% xlswrite('test',aczero_cell)
cd(outpath); 
ACdata_vs_model = merge_cells(ac0data_cell, aczero_cell, 2); 
ACdata_vs_model = merge_cells(ACdata_vs_model, ac1data_cell, 2); 
ACdata_vs_model = merge_cells(ACdata_vs_model, acone_cell, 2); 
if ~isunix
    xlswrite('ACdata_vs_model',ACdata_vs_model);
else
    savecell(ACdata_vs_model,[],outpath,'ACdata_vs_model');
end
% savecell(std_outcell_full,[],outpath,['Est Result', '_', mode_filename]);

cd(cucd);

%return
% make table in LaTeX
%%
znames_long = cell(length(znames)*2, 1); 
for ii = 1:length(znames); 
    znames_long{ii*2-1} = znames{ii}; 
    znames_long{ii*2} =  ' '; 
end
tex_only = 1; 
MakeLaTeX_table(ac_data(:, :, 1), znames, znames, 'Data: $\\rho_0$', ...
    'Data: Contemporaneous Cross Correlation', 'ac0data', tex_only, outpath);

MakeLaTeX_table(aczero_cell(3:end, 2:end), znames_long, znames, 'Model: $\\rho_0$', ...
    'Model: Contemporaneous Cross Correlation', 'ac0model', tex_only, outpath);


MakeLaTeX_table(ac_data(:, :, 2), znames, znames, 'Data: $\\rho_1$', ...
    'Data: First Order Cross Correlation', 'ac1data', tex_only, outpath);

MakeLaTeX_table(acone_cell(3:end, 2:end), znames_long, znames, 'Model: $\\rho_1$', ...
    'Model: First Order Cross Correlation', 'ac1model', tex_only, outpath);

cd(outpath);

fname = 'AC_DataVsModel';
f1=fopen([fname, '.tex'],'w');
         fprintf(f1, '\\documentclass{article} \n');
         fprintf(f1, '\\usepackage[cm]{fullpage} \n'); 
         fprintf(f1, '\\begin{document} \n');
         fprintf(f1, '\\input{ac0data.tex} \n');
         fprintf(f1, '\\vspace{3 mm} \n');
         
         fprintf(f1, '\\input{ac0model.tex} \n');
         fprintf(f1, '\\vspace{3 mm} \n');
         
         fprintf(f1, '\\input{ac1data.tex} \n');
         fprintf(f1, '\\vspace{3 mm} \n');
         
         fprintf(f1, '\\input{ac1model.tex} \n');
         fprintf(f1, '\\vspace{3 mm} \n');
         
        fprintf(f1, '\\end{document} \n'); 

fclose(f1);
if ~isunix
   status = system(['pdflatex ', fname,'.tex', ' -quiet']);
else
    status = system(['pdflatex "', fname,'.tex', '" -quiet']);
end
   delete([fname,'.aux'], [fname,'.log'])

cd(cucd); 

cd(outpath)
save ac_workspace ac_data aczero_cell acone_cell znames znames_long
cd(cucd)

%return
%%
if ~isempty(acsimz_mat);
    clear med_acsimz_mat;
    for ii = 1:ana_out.num_valid + 1;
        if ii == 1;
            med_acsimz_mat(:,:,:,ii)=ac_data;
        else
            med_acsimz_mat(:,:,:,ii)=median(acsimz_mat(:, :, :, :, ii-1),4);
        end
    end
    plotmultiple_acorr(med_acsimz_mat, in_ac);
end

%%
% 4.2) of the selected states
in_ac.modenames = ana_out.modelnames;
in_ac.varnames  = ana_in.state_sel;
in_ac.savepath = save_path;
in_ac.savename = 'Sel States AC Median';
 in_ac.flag_block = 1; 
if ~isempty(ana_out.acsims_mat)
    clear med_acsims_mat
    for ii = 1:ana_out.num_valid;
        med_acsims_mat(:,:,:,ii)=median(ana_out.acsims_mat(:, :, :, :, ii),4);
    end
    plotmultiple_acorr(med_acsims_mat, in_ac);
end
%%
% 4.3) full cross-correlogram vs. data
in_acfull.modenames =  [ana_out.modelnames(:); 'data'];

%varnames_correl={'u','gdp','C','I','v','find','sep'};
if exist('varnames_correl','var')==0
    in_acfull.varnames  = znames;
else
    in_acfull.varnames  =varnames_correl;
end
in_acfull.savepath  = save_path;
in_acfull.titfont   =12;
in_acfull.savename  =['Obs AC full ',subfname(2:end)];

%in_acfull.title = {'Obs AC full'};
if ~isempty(acsimz_mat);
    plotmultiple_acall(acsimz_mat, ac_data, in_acfull);
end